#!/usr/bin/python3.9
# -*- coding: utf-8 -*-
# @Time    : 2021/10/7 19:05
# @Author  : YHSimon

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import linear_model

'''
分析难点：
（1）hypothesis = X * theta.T
X = [x0 x1 x2 ... xn]，其中x0~xn是m*1的列向量，所有X是m*(n+1)的矩阵；
X是包含所有训练数据的特征矩阵，每一行都代表一条数据
theta = [theta0 theta1 theta2 ... theta-n]，1*(n+1)参数向量
所以hypothesis的值是m*1的向量，包含所有测试值的向量
（2）theta = theta - (alpha/y.size) * ((X * theta.T - y).T * X)
根据上面的分析X * theta.T是m*1的向量，包含所有测试值的向量；
(X * theta.T - y).T就是1*m的行向量，X是m*(n+1)的矩阵，故相乘结果为1*(n+1)的向量；
'''


def computeCost(X, y, theta):
    inner = np.power(((X * theta.T) - y), 2)
    return np.sum(inner) / (2 * len(X))


def gradientDescent(X, y, theta, alpha, epoch):
    """return theta,cost"""

    tmp = np.matrix(np.zeros(theta.shape))
    parameters = int(theta.flatten().shape[1])  # 参数θ的数量
    cost = np.zeros(epoch)  # 初始化一个ndarray，包含每次epoch的cost
    m = X.shape[0]  # 样本数量m

    for i in range(epoch):
        # 利用向量化一步求解
        tmp = theta - (alpha / m) * (X * theta.T - y).T * X

        # 以下是不用Vectorization求解梯度下降
        #         error = (X * theta.T) - y  # (97, 1)

        #         for j in range(parameters):
        #             term = np.multiply(error, X[:,j])  # (97, 1)
        #             temp[0,j] = theta[0,j] - ((alpha / m) * np.sum(term))  # (1,1)
        theta = tmp
        cost[i] = computeCost(X, y, theta)
    return theta, cost


def show():
    path = 'ex1data1.txt'
    # names添加列名，header用指定的行作为标题，若无标题且指定标题则设为None
    data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
    print(data.head())
    print(data.describe())
    print(data.get('Population'))
    # plt.figure(figsize=(8, 5))
    plt.scatter(data.get('Population'), data.get('Profit'))
    # data.plot(kind='scatter', x='Population', y='Profit', figsize=(8, 5))
    plt.show()


if __name__ == '__main__':
    path = 'ex1data1.txt'
    # names添加列名，header用指定的行作为标题，若无标题且指定标题则设为None
    data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
    data.insert(0, 'Ones', 1)
    cols = data.shape[1]  # 列数
    X = data.iloc[:, 0:cols - 1]
    y = data.iloc[:, cols - 1:cols]  # 取最后一列，即目标向量
    print(X.head())
    print(y.head())

    X = np.matrix(X.values)
    y = np.matrix(y.values)
    theta = np.matrix([0, 0])  # theta是一个(1,2)矩阵
    print(X.shape, y.shape, theta.shape)
    alpha = 0.01
    epoch = 1000
    final_theta, cost = gradientDescent(X, y, theta, alpha, epoch)
    # 使用拟合的参数计算训练模型的代价函数(误差)
    computeCost(X, y, final_theta)

    # 绘制线性模型以及数据
    # np.linspace()在指定的间隔内返回均匀间隔的数字
    x = np.linspace(data.Population.min(), data.Population.max(), 100)  # 横坐标
    f = final_theta[0, 0] + (final_theta[0, 1] * x)

    fig, ax = plt.subplots(figsize=(6, 4))
    ax.plot(x, f, 'r', label='Prediction')
    ax.scatter(data['Population'], data.Profit, label="Training Data")
    ax.legend(loc=2)
    ax.set_xlabel('Population')
    ax.set_ylabel('Profit')
    ax.set_title('Predict Profit vs. Population Size')
    plt.show()

    # 绘制代价函数与迭代次数
    fig1, ax = plt.subplots(figsize=(8, 4))
    ax.plot(np.arange(epoch), cost, 'r')
    ax.set_xlabel('Iterations')
    ax.set_ylabel('Error vs. Training Epoch')
    plt.savefig('代价与迭代次数的关系.png')
    plt.show()

    # (3) 预测人口规模为3.5 和7万的利润值
    print(final_theta[0, 0], final_theta[0, 1])
    print((final_theta[0, 0] * 1 + (final_theta[0, 1] * 7)) * 10000)
    print('------------')

    # 使用scikit-learn的线性回归函数
    model = linear_model.LinearRegression()
    model.fit(X, y)
    x3 = np.array(X[:, 1].A1)
    f3 = model.predict(X).flatten()

    fig, ax = plt.subplots(figsize=(8, 5))
    ax.plot(x3, f3, 'r', label='Prediction')
    ax.scatter(data.Population, data.Profit, label='Training Data')
    ax.legend(loc=2)
    ax.set_xlabel('Population')
    ax.set_ylabel('Profit')
    ax.set_title('Predicted Profit vs. Population Size')
    plt.show()
